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Abstract 

The discrete time path integral Monte Carlo (PIMC) with a one-particle density matrix approx- 
imation is applied to study the quantum phase transition in the coupled double-well chain. To 
improve the convergence properties, the exact action for a single particle in a double well potential 
is used to construct the many-particle action. The algorithm is applied to the interacting quantum 
double- well chain for which the zero-temperature phase diagram is determined. The quantum phase 
transition is studied via finite-size scaling and the critical exponents are shown to be compatible 
with the classical two-dimensional (2D) Ising universality class - not only in the order-disorder 
limit (deep potential wells) but also in the displacive regime (shallow potential wells). 
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I. INTRODUCTION 



The two-level tunneling model provides a phenomenological description of the low- 
temperature properties of glassy materials [l|. In the simplest case, an isolated tunneling 
system can be represented by a particle moving in a double-well potential. Experimental 
findings have suggested that the interactions between the tunneling systems play a cru- 
cial role in the low temperature behavior which deviates from the predictions of the non- 
interacting two-level systems [2|. The model Hamiltonian of the system with L particles is 
then given by 




where X{ is the (one-dimensional) displacement of the z-th particle (i = 1, . . . , L) of mass 
/i from a reference position, = denotes the momentum operator, U(xi) is a local 
potential for the displacement of the i-th particle that is usually assumed to be a double well 
potential, and V(xi,Xj) describes the interaction between particles, see Fig. [IJ Apart from 
glassy materials, this coupled double-well model has been applied to other systems, including 
structural phase transitions of a wide range of systems, e.g. uniaxial ferroelectricsfl. Most 
numerical computations devoted to understanding the interacting double-well model have 
mainly treated the problem in the framework of the classical 4 model or have been limited 
in the "two-state" limit by studying the corresponding Ising model. These simplifications 
reveal the difficulties inherent in simulations of the quantum coupled double-well model. 

In this paper we present an efficient path integral Monte Carlo (PIMC) algorithm to study 
interacting particles, each of which is confined to a double well potential. The method is 
presented in the next section, and it is applied to the one-dimensional interacting double- well 
model in the Sec. IHH which also contains the results: the phase diagram and the discussion 
of the universality class of the quantum phase transition. 



II. THE METHOD 



The partition function of ([I]) is given by 

Z = [ dxp(x,x;/5) , (2) 
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where x = (xi, . . . , xl) is the displacement configuration of the whole system and 

p(x,x';/?) = (x|e-^|x') (3) 

is the density matrix, with j3 = 1/T the inverse temperature. Observables that are diagonal 
in the displacements, like the the m-th moment (x™), are given by 

(0(x)> = i J dxp(x,x;/3)0(x). (4) 

By splitting the Hamiltonian flTJ into two (non-commuting) parts H = Ha + Hb and using 
the Suzuki- Trotter identity, one arrives at the path integral formula for p: 

M 



e -f3(H A +H B ) _ J im 

Af^oo . 



(5) 



with r = (3/M. The conventional choice for Ha and Hb is the kinetic energy for 7"Ci (which 
is diagonal in the momenta and the potential plus interaction energy U+V for Hb (which 
is diagonal in the displacement variables). In the lowest-order of the commutator expansion 
- the so-called primitive approximation, the high-temperature density matrix becomes 

p(x, x'; r) = f[ e-W^+v^Mpaixi, x[; T ) e -^ u ^ +v ^) (6) 

where po(xi,x' i ;r) is the free particle density matrix. This choice leads to bad convergence 
properties in the Trotter number M {4] because of the fractal character of a trajectory of a 
free quantum mechanical particle described by the term Ha- 

The purpose of this paper is to demonstrate the efficiency of another choice for Ha 
and Hb by treating the single particle diffusion within a double well potential exactly and 
separately from the particle interactions. Doing this, we have 

l 2 

n A = £ ^ + u( Xi ) 
1=1 

H B = ^Vfaxj) . (7) 

i<3 

This strategy is expected to be most promising in the case when the interactions are much 
weaker than the mean potential energy of the particle. 

With (J7|) the path integral expression for the density matrix becomes 

„ M-l 

p(x,x;/5)= lim / dx 1 . . . dx M -i TT Pa(x™, x m+i ; r)p B (x m , x m+1 ; r) (8) 

J m=0 



with x = x = xm and 



L 



p A (x,x';r) 



i=l 



Pb(x,x';t) 




(9) 



i<j 



where 



P 



{1) (x,x';t) 



{x\e 2 »\x)e 



U(x) + U(x') 
1 o 
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(10) 



is the one-particle density matrix for a single particle in a potential U. For a double well 
potential this is not known analytically but can easily be computed numerically with the 
matrix multiplication method^]. This method is based on the recursion formula 



and the fact that in the limit St — > the one-particle density matrix p^ can be factorized 
into the kinetic and potential energy part: 



By squaring the density matrix k times, we will lower the temperature by a factor of 2 k 
and reach the required temperature r. For a given potential U(x), the limits of integration, 
x m j n and x max , are chosen appropriately - not too large for computational reasons and not 
too small for numerical accuracy. Once the limits are set, a fine grid between x m i n and 
^max should be constructed for the numerical integrations. The spacing between successive 
grid points should be sufficiently small to ensure the high accuracy. We store this one- 
particle density matrix in a two-dimensional array as a look-up table for use during the 
simulations, and employ a simple bilinear interpolation to determine the matrix elements 
for any point (x^x^) within [x min ,x max ] in the continuous position space. We note that the 
symmetric break up of the propagator in the form of Eq. (|T0|) satisfies a unitarity condition 
pW(5r)pW(— 5t) = 1, which can be utilized to reduce errors resulting from discretization of 
r, as discussed in jfj]. 

Path integral Monte Carlo means the evaluation of the integral (jHJ) via importance sam- 
pling of the configurations (xi, . . . ,xm-i) (for fixed M) with the appropriate weight given 
by the integrand of (jHJ). Here we use a single step update scheme: Let X = (x 1; . . . , ~Km-x) — 
{xi tm } be the current configuration. We generate a new configuration X' which differs from 




(11) 




(12) 
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the old configuration X only in a single particle displacement in a particular time slice: 
x' im = Xi t7n + 6, where 6 G [— e, e] is a uniformly distributed random number, and e the step 
size. The acceptance probability w(K — > X') of this new configuration should be chosen to 
fulfill detailed balance with respect to the weights of of the old and new configurations; a 
possible choice is 



w(X -> X) 



mm 



mm 



M-l 



j TT Pi(Vli x mi 1 ')PB( x m ! X m+l! T ) 



m=0 



PA(x m -i,x m ; r)p B (x m ,x m+ i; r) 



^ g — tAV(X,X') \ X i,rn-li x i : rn' T )P^ ^ ( X i,m' X i,m+li r ) 

P i, x i,m— 1) x i,m] T)p^' (li.mj 2-i,m+l) 



(13) 



III. THE ONE-DIMENSIONAL MODEL AND RESULTS 



To test the above algorithm we focus here on a one- dimensional geometry in which par- 
ticles interact only with their nearest neighbors, and assume V(xi,Xi+i) to be quadratic in 
the displacements 

L 

Kb = ^2 V ( Xi > x ^ = ~ JiXiXi +^ ' ( 14 ) 

i<j i=l 

(cf. Fig. [2]). Furthermore we choose homogeneous interaction strength Jj = J > 0, and the 
double- well potential in the symmetrical form with two minima located at ±1: 

U(x) = V Q (x A -2x 2 ). (15) 

Periodic boundary conditions are imposed. 

The model (PQ) with (1141) and (I15p has a Z 2 -symmetry (x, — > — XjVz) and corresponds to 
a quantum version of a 4 theory, which is expected to belong to the universality class of 
1 + 1-dimensional Ising model. Suppose that the height of the potential barrier Vq is large 
compared to the energy scale of the particle executing small oscillations in one of the double 
wells. The model is then equivalent to the one-dimensional Ising model in a transverse field 

H TIM = -r5>f- J^Vftrj, (16) 



where the transverse field T corresponds to the tunneling splitting in the double-well pro 
lem. Therefore we expect it to display a zero-temperature quantum phase transition 
from a disordered phase with (xj) = to an ordered phase with (x,) ^ at a critical inter- 
action strength J c (for fixed p and Vq). According to the universality hypothesis, the same 



universality class, i.e. the 2D Ising, should extent to the region where Vp is small compared 
to the interactions between particles, the so-called displacive region 3, 7, sj]. 

For a given value of the parameters Vq and J we computed the following quantities: the 
average of the displacement m (i.e. the magnetization in the spin formulation), defined as 
m(L, M) = J2i 12n ( x i,n) , where x itn is the position of the i-th particle at the time step 
n with respect to the zero position of the local potential Vd w ; the fourth-order cumulant of 
the magnetization given by g = | ^3 — (m 4 ) / (m 2 ) 2 ^, where (- • ■ ) denotes the expectation 
value over MC configurations; the susceptibility defined as x — Lf3 ((m 2 ) — (|m| 2 )). Close 
to a quantum critical point, one expects [Qf] observables O to scale as 

O = L~ x °0 (d L 1 ^, y-\ (17) 



where xo is the scaling dimension of the observable O, v the correlation length exponent 
and z the dynamical exponent. If the transition falls into the Ising universality class, the 
dynamical exponent z is unity [9]. In the following we assume this to be the case and check 
whether our data are compatible with this. We choose a fixed value of the aspect ratio L//3, 
corresponding to z = 1, so that the finite-size scaling function of these quantities involves 
only one variable, i.e. O = L- x °0(5 L l / U ). For a given Vo, the deviation from the critical 
point is parameterized as 5 = J — J c . The scaling dimension is given by x m = —f3 m /v 
for the magnetization \m\, x x = 7/1/ for the magnetic susceptibility x an d x g = for the 
dimensionless fourth cumulant g. Typically we executed 10 6 — 10 7 MC steps to thermalize the 
system. Once in equilibrium, we generated 4 — 5 x 10 7 MC configurations for measurements, 
which were carried out every 5 MC steps. We have considered a wide rage of values of Vq 
between 0.01 and 5. At a fixed value of Vo we varied the strength of the ferromagnetic 
interaction J for system sizes up to L = 64 to localize the zero-temperature critical point 
J c and to carry out the finite-size analysis. 

To confirm the accuracy of the one-particle density matrices calculated by matrix mul- 
tiplication method, we first compare the distribution of displacements of the particles in 
the absence of the interaction obtained by PIMC with the distribution calculated by solv- 
ing numerically the single particle Schrodinger equations. For the latter, we calculated the 
first N = 50 energy eigenvalues E n and the corresponding eigenstates ip n (x)', the distribution 
function of the displacement is then computed by P{x) = Xm=i l ? /'™( x )| 2 e~ fiEn j ^2 n=1 e~ fiE " n . 
As shown in Fig. [3] for (3 = 16 by using r = 0.25, the excellent agreement confirms the ac- 
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curacy of the density matrices. Furthermore, we compare the results from PIMC within the 
one-particle density matrix approximation with those in the primitive approximation for the 
same parameters (e.g. Vq and J), as shown in Fig [5] as a typical example for the depen- 
dence of magnetization \m\ and its fourth-order cumulant g on the time step r. The results 
presented are averaged over 16 samples for each time step. We find that, with r = 0.05, 
the results obtained by the primitive approximation converge to the same values computed 
with the one-particle density matrix for r < 0.25 within the statistical error bars. The CPU 
time required on an Intel Pentium processor (2.40GHz) to calculate 500000 MC steps for a 
system size L = 32 and (3 — 16 within the one-particle density matrix approximation using 
t = 0.25 is about 780 seconds, and with the primitive approximation using r = 0.05 is 
about 1485 seconds. The efficiency of the calculations with the one-particle density matrix 
is gained from the fast convergence with respect to the number of time slice. After a careful 
examination, we are convinced that the time step r = 0.25, used in our simulations for 
the high-temperature one-particle density matrix, is sufficiently small for the convergence. 
We carried out 8 iterations for the matrix multiplication to generate a one-particle density 
matrix with r = 0.25 and the spacing between neighboring points within [x m i n , x max ] was 
chosen to be 0.01. In all cases studied we used a wide interval of [x min , x max ], e.g. [—10, +10] 
for Vq — 3, for the iterative integrations and then truncated this interval to a smaller one 
while storing into the look-up table for PIMC simulations. The appropriate values for the 
interval [x m i n , x max ] in the look-up table were justified by doing a short run of the PIMC 
simulation to check whether the particles would move beyond the chosen boundaries. 

In Fig. [5] we present the zero-temperature phase diagram, in which the critical value J c 
is estimated by the intersection of g(J) curves at a given V for various system sizes (up to 
L = 64) with fixed aspect ratio L//3 = 2. We note the lack of monotonicity of the critical 
J c with respect to the potential barrier Vo; J c decreases with V in the deep well region, 
while it increases with Vo in the small Vo region. This non-monotonic behavior of J c {Vo) is 
qualitatively reproduced within the mean field approximation: Consider the effective single- 
site Hamiltonian including a mean-field term 

2 

H m{ = ^- + V (x 4 - 2x 2 ) - 2Jxm, (18) 
where the order parameter m is the expectation value of the displacement x in the ground 
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state if)o(x, m) of 7Y m f and is determined self-consistently via 

m = dxx\ipo(x,m)\ 2 . (19) 



Varying J and solving the non-linear equation (fl9|) for m numerically, the critical point can 
be estimated as the value of J above which a non-zero solution exists. 

The mean-field result for J c , depicted in Fig. [5j shows the same non-monotonic behavior 
as our results for J c from the PIMC and has a maximum at Vq = 1. This behavior of J c 
can be understood as follows: In the region Vq ^> 1 the potential has two deep minima 
separated by a barrier Vq giving rise to a nearly degenerated ground state doublet that is 
well separated from the rest of the spectrum, as shown in the upper panel of Fig. [51 This 
is called the order-disorder limit [3), [h]] in the interacting double- well model. The energy 
difference between the ground state and the first excited state, i.e. the tunneling splitting 
T, is reduced as Vq grows, which results in a decrease of the critical ordering term J c . In 
the region Vq <C 1, on the other hand, the potential has two shallow minima and the two 
lowest energy levels are not well separated from the rest of the spectrum. This is called the 



displacive regime [111] in the interacting double-well model. In this displacive region, the zero 
point energy of a single particle lies above the barrier of the local potential so that the local 
potential is effectively in a single-well form. Without the interaction term, the particles 
fluctuate around the x = position (cf. Fig. [3]); switching on the interaction shifts the 
displacement expectation value (x) away from the origin and at the critical coupling J c the 
systems undergoes a displacive transition from a symmetric (disordered) phase to a broken 
symmetry (ordered) phase. The key factor for the strength of the critical displacing force 
J c in this case is the width of the local potential, which decreases with increasing Vo: the 
wider the local potential, the weaker the force J needed for the displacement. Therefore, 
the critical value J c increases with Vq in the displacive regime. 

For a particular value of Vq we can use the scaling form given in Eq. (ITT)) for g, \m\ and 
X to extract values of the critical exponents. In all cases a good data collapse is achieved 
with the exponents (3 m = 1/8, 7 = 7/4 and v = 1.0, which is representative of the classical 
2D Ising universality class. In Fig. [6] we show the finite-size scaling plots for Vq = 3 and 
Vq = 0.01. For Vq = 0.01 which is well inside the displacive regime, the quality of the scaling 
decreases and corrections to scaling become more pronounced. Interestingly, the peak of the 
scaling function x(t) of the susceptibility is shifted away from t — for Vq — 0.01, whereas 
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it is at t = for Vq = 3, indicating the non-universality of the scaling function. Our results 
for Vo 1 indicate that the model is in the Ising universality class even in the displacive 
regime. For the interaction in the form given in Eq. (114)) . we expect that a phase transition 
in the same universality still occurs when the local potential is reduced to only a quartic 
term. To provide support for this we carried out simulations for the model with a local 
potential given by U (x) = x 4 which exhibits a single well. Our results depicted in Fig. [7] 
suggest that the critical behavior of this one-well model is indeed consistent with 2D Ising 



universality 



121 ] . We note that the coupling term (fll]) that we use can be brought into a 



form that is more reminiscent of a lattice version of the standard 4 (quantum) field theory: 

L L 

2 



Kb = J - - ^+i) 2 - J x * ' ( 2 °) 



i=i i=i 



Together with the local potential ( fl5l) this implies that the corresponding continuum model 
for a scalar field <p contains a (d<p/dx) 2 -term and a potential energy of the form Vo[0 4 — 
(2 + J/Vq)(/) 2 )]. Since J is always positive and can be made arbitrarily large, this model 
has always a phase transition (at zero temperature). On the other hand, the field theory 
with a potential energy that has only a single minimum, like the pure quartic potential 
Vo0 4 , does not have a phase transition. We checked, within mean-field as well as with PIMC 
simulations, that the corresponding lattice model 

i=i ^ * 

also does not have a phase transition. 

To summarize, we have demonstrated that PIMC within the one-particle density matrix 
approximation is an efficient method to simulate quantum interacting many-body systems, 
in which particles are confined in a local potential and interact with each other. Using this 
method we have studied the zero-temperature phase transition of the coupled double-well 
chain, both in the order-disorder case, corresponding to a coupled two-level tunneling system, 
and in the displacive regime, in which the interaction dominates over the double-well struc- 
ture. Based on this numerical scheme, our further study will include the double-well model 
coupled through long-range/random interactions and coupled to a dissipative bathjis), 
In the presence of quenched disorder in the coupling, even for the case without dissipation, 
implementation of many improved PIMC methods, e.g. Fourier PIMC techniques or cluster 
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algorithms, becomes complex and the computational efficiency reduces. This motivates the 
choice of a method which provides easy performance and can be extended to the random 
case in a straightforward way, as the technique applied in this paper does. 
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FIG. 2: Representation of the model denned in Eq. (pQ) in the one-dimensional form. 
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FIG. 3: The distribution of displacements of the particles, calculated for L = 32 and j3 = 16, 
with J = 0. A comparison with numerical solutions (indicated by the solid line) of the Schrodinger 
equation is shown. The excellent agreement confirms the accuracy of the density matrices. 
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FIG. 4: The magnetization (\m\) and its fourth-order cumulant g, calculated by using the one- 
particle density matrix and the primitive approximation. The model parameters are Vo = 1 and 
J = 0.56 for a system size L = 32 at temperature (3 = 16. The values obtained from both methods 
are compatible in the small r limit. 
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FIG. 5: Lower panel: The phase diagram of the coupled double-well chain: the critical ratio 
J c /Vb as well as the critical interaction J c as functions of the depth of the potential well Vb; the 
ordered phase is located above the curves and the disordered phase is below the curves. The critical 
J c obtained by the mean field approach is indicated by the dashed line. Upper panel: The low 
lying energy eigenvalues of the one-particle Hamiltonian for various Vo, determined by numerical 
solutions of the Schrodinger equation. The dotted line at E = indicates the top of the potential 
barrier. 
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FIG. 6: Scaling plots of the cumulant(a), the magnetization(b) and the susceptibility(c) for 
Vq = 0.01 (left) and Vq = 3 (right) using the two-dimensional Ising universality. 
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FIG. 7: Scaling plots of the cumulant(a), the magnetization (b) and the susceptibility(c) for the 
single- well model with U (x) = x 4 using the two-dimensional Ising universality. In the scaling plot 
(c), the system sizes only range from L = 16 to L = 64. 
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